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Abstract 

We address the relationship between membrane microheterogeneity and anomalous subdif- 
fusion in cell membranes by carrying out Monte Carlo simulations of two-component lipid 
membranes. We find that near-critical fluctuations in the membrane lead to transient subdif- 
fusion, while membrane-cytoskeleton interaction strongly affects phase separation, enhances 
subdiffusion, and eventually leads to hop diffusion of lipids. Thus, we present a minimum 
realistic model for membrane rafts showing the features of both microscopic phase separation 
and subdiffusion. 

PACS: 87.14.Cc, 87.15.ak, 87.16.dj, 87.16.dt. 

Keywords: Lipid membranes, DMPC, DSPC, phase separation, membrane rafts, subdiffu- 
sion. 
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Introduction 



Anomalous subdiffusion in cell membranes is an intriguing phenomenon whose molecular ori- 
gins are still a subject of debate [TH3]. From the general viewpoint, the phenomenon is related 
to (dynamic) microheterogeneities in the properties of the cell membrane creating a rugged 
energy landscape for the diffusing protein or lipid molecule. The current understanding 
of membrane microdomains embraces the concepts of lipid rafts [4] and cytoskeleton-based 
picket fence [5j . The general concept of membrane rafts [6] implies a dynamic interplay be- 
tween the membrane local composition and phase, and local membrane-protein interactions, 
which can result in local lipid demixing and phase separation. This view is supported by the 
recent experimental evidence demonstrating that local lipid demixing and phase separation 
can be induced by crosslinking [7], change in the local membrane curvature [8]j, and inter- 
action with a cytoskeleton [9j. Additionally, the presence of the membrane-associated actin 
network leads to a shift in the membrane phase transition temperature and, potentially, to 
broadening of the phase transition [9]. 

The efficiency of external perturbations in changing the local properties of the mem- 
brane should be strongly enhanced in the vicinity of the membrane critical point. This is 
indeed supported by recent experimental observations [10J. Moreover, it was suggested that 
dynamic microdomains in the cell membrane are nothing but near-critical fluctuations in 
the local composition and phase of the membrane [TlJ. All this is in agreement with the 
observation that the composition of the cell membrane is adjusted to keep it above the phase 
transition temperature, and is regulated in response to environmental changes to maintain 
this condition [T2] . In cases where this mechanism fails, cold shock damage takes place [T3] . 

In this paper, we address the relationship between membrane microheterogeneity and 
anomalous subdiffusion in cell membranes by carrying out Monte Carlo simulations of two- 
component lipid membranes on experimentally relevant spatial scales (~1 /im) and time 
intervals (~1 s), with a special emphasis on interactions with the model membrane skeleton. 

We demonstrate that (i) near-critical fluctuations in a free lipid membrane can lead to 
transient anomalous subdiffusion, (ii) phase separation in two-component (and, hence, mul- 
ticomponent) lipid membranes can be strongly affected by interaction with the membrane 
skeleton, which, depending on the temperature and membrane composition, can either lead 
to precipitation of highly dynamic membrane domains (rafts), or prevent large-scale phase 
separation, and (iii) interaction with the membrane skeleton enhances anomalous subdiffu- 
sion and eventually leads to hop-diffusion of lipids. Thus, we construct a minimum realistic 
model for membrane rafts showing the features of both microscopic phase separation and 
anomalous subdiffusion. 

The binary lipid system consisting of l,2-dimyristoyl-sn-glycero-3-phosphocholine 
(DMPC) and l,2-distearoyl-sn-glycero-3-phosphocholine (DSPC) is chosen as a generic 
model of a two-component membrane. The particular choice of this system is motivated 
by the fact that it is well studied both in vitro fTH - [2T] and in silico (see [20H22] and refs. 
therein) . Even though the exact temperature range of phase separation in this membrane is 
too high for most living organisms, the main qualitative conclusions of the study related to 
behavior of cell membranes are not affected by the particular thermodynamic parameters of 
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the lipid mixture. 

From the physical perspective, it is assumed that living organisms try to maintain the 
composition of cell membranes to keep them above the phase transition temperature [T2] . 
Therefore, to be biologically relevant, in the present paper we mostly focus on the upper 
part of the phase diagram. 

Materials and Methods 

Experimental 

The saturated phospholipids DMPC and DSPC with the melting temperatures 297 and 328 
K were purchased from Avanti Polar Lipids (Alabaster, AL). Multilamellar vesicle (MLV) 
suspensions were obtained by rehydration of a dry lipid film. Excess heat capacity curves 
of DMPC/DSPC MLV suspensions in a 10 mM Hepes buffer, pH 7.4, were obtained using a 
VP-DSC calorimeter (MicroCal, Northampton, MA) at a scan rate of 2 3 K/h. 

The empirical phase diagram based on the temperatures of the onset and completion of 
the phase transition determined from the scanning calorimetry data as described [20l [21] are 
in agreement with previously reported results (see Supporting Material). 

Monte Carlo simulations 

Our approach to lattice-based Monte Carlo (MC) simulations of a two-component membrane 
is generally similar to the one described previously [20l [21]. To facilitate efficient simula- 
tions on experimentally relevant spatial scales (~1 fim) and time intervals (~1 s), where a 
particular type of lipid packing and fine molecular details should be of little importance, we 
further simplified the model and represented the membrane as a square lattice, each node 
of which represents a molecule of one of the two lipid types, which can be in either gel or 
fluid conformational state. An elementary MC step consists of two substeps: (i) an attempt 
to change the state of a randomly chosen lipid and (ii) an attempt of position exchange of a 
randomly chosen next-neighbor pair of lipids. As in [21], a rate function was introduced for 
the next-neighbor exchange step to ensure ~40 times slower lipid diffusion in the gel phase. 
For anLxL lattice, an MC cycle consists of a chain of L 2 elementary MC steps. For every 
lipid composition and temperature studied, the membrane was first equilibrated, if required, 
and simulations were run for (6 — 20) x 10 6 MC cycles to collect the necessary data. More 
details on the simulation procedure are given in the Supporting Material. 

Simulations were carried out on an L x L = 600 x 600 (400 x 400 when modeling the 
effects of the membrane skeleton) square lattice with periodic boundary conditions. By 
assuming the average lipid headgroup size of 0.8 nm, this corresponds to a membrane with 
an experimentally relevant size of 0.48 x 0.48 /xm 2 (0.32 x 0.32 /xm 2 ). By comparing the 
DMPC diffusion coefficient obtained in our simulations of pure DMPC in the fluid phase 
at 304 K with its experimental values at the same temperature (3 — 6 /xm 2 /s f23H25] ). we 
found that one MC cycle corresponds to ^50 ns, and thus, our simulations cover processes 
on timescales up to ^0.3-1 s. 
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Model for the membrane skeleton 

The membrane skeleton was modeled by a random Voronoi tessellation satisfying the pe- 
riodic boundary conditions. For simulations with L = 400, random tessellations with 
N = 36 compartments were used, which gives the average linear size of the compartment 
£ = LN~ 1/2 = 66.6 lattice units « 53 nm. The generated filament meshwork was projected 
onto the square lattice thereby creating a set of pixels representing the locations of the fila- 
ments. Each of these locations can be assigned to be a cytoskeleton pinning site. To mimic 
the membrane-cortical skeleton interaction, a simple rule inspired by experimental data on 
lipid interactions with transmembrane proteins [26] was followed: a lipid located at a fila- 
ment pinning site is forced to assume the gel conformation with no explicit restrictions on its 
mobility. Thus, the pinning site does not present an obstacle for lipid diffusion: its effect on 
diffusion is indirect and takes place solely due to a lower lipid mobility in the gel-state local 
environment. The effect of the varying strength of the membrane-skeleton interaction was 
modeled by randomly assigning a fraction of filament position pixels to be filament pinning 
sites. Simulations were carried out with the filament pinning density set to 25%, 50%, and 
100%; in these cases the total number of pinning sites amounted to approximately 1%, 2%, 
and 4% of the total membrane area. 

Our choice of immobile pinning sites is justified by experimental observations demon- 
strating that band 3 and ankyrin strongly bound to the membrane skeleton show a very low 
diffusion coefficient of ~ 10 -4 — 10 -3 //m 2 /s over time intervals up to tens of seconds [271(28] . 

It should be pointed out that the approach used in the present work to account for 
interactions of lipid molecules with membrane proteins is conceptually similar, though not 
identical, to the one used in several previous MC simulation-based studies [29H32] . In these 
works, it was assumed that proteins, preferentially wetted by one of the membrane phases, 
could freely diffuse in the membrane, which results in their accumulation in this phase and, 
in case of two-component membranes, in fact enhances large-scale phase separation. (The 
situation is more complicated in case of active proteins - for details, see [32] and refs. therein.) 
It was found out that, under certain conditions, interaction of a single-component membrane 
with small transmembrane proteins or peptides can even lead to emergence of a closed loop 
of gel-fluid coexistence with a lower critical point [33J. It should be pointed out, however, 
that these works focused only on structural properties of protein-loaded membranes, and did 
not address the issues of diffusion in the membrane. 

What sets the present work apart from the above-mentioned studies, is that here we 
consider interaction of the membrane with the immobile cytoskeleton and study its effects 
on phase separation and diffusion of lipids in the membrane. 

Analysis of lipid diffusion data 

Positions of a small fraction of lipid molecules (150 for simulations with L = 600 and 50 for 
simulations with L = 400, which amounted to 0.04% and 0.03%, respectively) were recorded, 
and the time- and ensemble-averaged mean-square displacement (MSD) was determined as 
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follows: 



MSD(r) = — — ^ {[*(*) ~ x(t + r)f + [y(t) - y(t + r)] 2 } , (1) 

where t, i, and £ max are times measured in units of MC cycles; here, r denotes the time lag, 
and i max is the total length of the lipid molecule trajectory. No difference between time- and 
ensemble-averaged MSD (Eq. [T]) and time-only- and ensemble-only-averaged MSDs beyond 
the statistical error level was observed. 

In case of normal diffusion, the MSD grows in a linear fashion with time: MSD(r) = ADr 
in 2D, where D is the translational diffusion coefficient. In case of anomalous subdiffusion, 
the MSD shows a slower sublinear power-law growth MSD{r) ~ with < /3 < 1 (see, e.g., 
[34]). An alternative description of diffusion showing deviations from the normal behavior 
(also in cases where it cannot be described in terms of subdiffusion) can be provided using 
an effective time-dependent diffusion coefficient D{r). 

Therefore, to characterize the behavior of the MSD curves, the local exponent of the 
mean-square displacement 

/3msd(t) = d\ogMSD(r)/d\ogT, (2) 
and the effective time-dependent diffusion coefficient 

D(t) = htMSD{r)/dr (3) 
were calculated. 



Simulation of FCS experiments 

To simulate fluorescence correlation spectroscopy (FCS) [35] measurements, the tracked par- 
ticles (as above, in the amount of 150 for simulations with L = 600 and 50 for simulations 
with L = 400) were assumed to be fluorescent. Fluorescence intensity fluctuations 5F(t) = 
F{t) — (F) about the mean intensity (F) in a 2D Gaussian detection spot exp(— 2r 2 /r$) were 
recorded, and their autocorrelation function G(r) = (SF(t)5F(t + r))/(F) 2 was calculated. 
The detection spot size was set to r = 31 lattice units « 25 nm, the size experimentally 
achievable using the stimulated emission depletion (STED) FCS technique [36j. We addi- 
tionally note that this detection spot is much smaller than the lattice size (vq/L ~ 0.05 for 
L = 600 and r /L w 0.08 for L = 400), which allowed us to avoid artifacts in the fluorescence 
autocorrelation function and additionally ensured that our simulations are experimentally 
relevant. 

FCS curves were averaged over nine different positions on the lattice; when studying the 
effects of the membrane skeleton, the results were additionally averaged over five random 
realizations of the filament network. 

Simulated FCS curves were analyzed using the model 

G(r) = G(0)/ [1 + (r/Tb)*"] . (4) 
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For /3 F cs = 1 this expression corresponds to normal diffusion, while for < /3pcs < 1 it 
provides a simple way to describe anomalous subdiffusion in FCS [35]. In case of normal 
diffusion, td is related to the diffusion coefficient of fluorescent particles D and detection 
spot size r : td = Tq/(AD). 

Notice that since FCS is sensitive not only to the MSD(t), but also to the higher moments 
of the distribution of displacements of a diffusing particle, the connection between G(r) and 
MSD(r) is straightforward only in case of a Gaussian distribution of displacements [35] . As 
a result, in case of anomalous diffusion of particles, generally /3pcs 7^ /5msd- 

Results and Discussion 

Phase and component separation in the membrane 

With appropriate tuning of lipid interaction parameters, the empirical heat capacity-based 
phase diagram obtained from our MC simulation is in agreement with our experimental data 
(Fig. [I]), as well as with previously published experimental and simulation data on the same 
lipid system (see Supporting Material). At the same time, the empirical phase diagram 
constructed on the basis of heat capacity data may differ from the real phase diagram of the 
system and thus not provide an insight into the microscopic structure and the dynamics of 
the membrane. 

More details on the phase diagram can be obtained based on binodal and spinodal curves 
of the system. The binodal curves, also known as coexistence curves, are the boundaries of the 
region in the phase diagram in which the equilibrated system shows a complete separation of 
the two phases. The spinodal encloses the region where the mixture is unstable with respect 
to small local fluctuations of the composition and always lies inside the area enclosed by the 
binodal, with the exception of a critical point, where the binodal and spinodal touch. 

To reconstruct the binodals and spinodals of the DMPC/DSPC lipid mixture, we ana- 
lyzed the static structure factors [37J for lipids Si,(k) and lipid states Ss(k) (see Supporting 
Material), reflecting the character of spatial fluctuations of the membrane composition and 
state. 

We found that outside the phase coexistence region, the Ornstein-Zernike (OZ) approx- 
imation [37J S oz (k) = 5(0)/ [1 + (&£) 2 ], where k is a wavenumber, and £ is a correlation 
length, provided an excellent description of the lipid-state structure factors Ss(k). The tem- 
perature dependences of £s(0) and £s were used to estimate spinodal and binodal curves. 
In particular, the spinodals were determined by extrapolating l/£s(0) dependence to zero 
crossing [371 [38], and binodals were determined from the condition d(l/£g)/dX = [38] . 
Interestingly, in some regions of the phase diagram, the binodal deviates quite strongly from 
the phase coexistence boundary estimated from the excess heat capacity data (Fig. [T]). This 
discrepancy indeed shows that the analysis of heat capacity data does not necessarily recover 
the real phase diagram of the system. Here, this behavior reflects the continuous charac- 
ter of the phase transition in the membrane, whose proper description requires a combined 
component and state phase diagram [39] . 

What process is reflected by the excess heat capacity curves where they fail to describe 
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Figure 1: Left-hand panel: Component and lipid state separation phase diagram of 
DMPC/DSPC lipid mixtures. Phase transition temperatures as determined from experi- 
mental (A) and simulated (o) excess heat capacity curves. Lipid state spinodal ( ), lipid 

state binodal ( ), and lipid demixing curves ( — ). Right-hand panel: representative snap- 
shots of equilibrium membrane configurations at temperatures and membrane compositions 
corresponding to filled squares in the left-hand panel. Lattice size: 600 x 600; scale bar: 200 
lattice units « 160 nm. 



the fluid-gel phase separation? An analysis of the structure factors Si,(k) for the lipid species 
helps to answer this question. It appears that, generally, two OZ components are required 
to describe these data: 5l(&) = S®±(k) + S^ify. Parameters of component 1 only weakly 
depend on the temperature and reflect demixing of lipids in the same state. Parameters 
of the second component Sl2(0) and £l2 show strong temperature dependences similar to 
those of £s(0) and £s, and thus describe appearance of dynamic microscopic domains, which 
cannot be treated as distinct thermodynamically stable phases [22J. We therefore define 
the lipid demixing curves as a set of points on the phase diagram satisfying the condition 
5li(0) = Sl2(0). The fact that the lipid demixing curves are in excellent agreement with 
the calorimetry-based empirical phase diagrams, supports the above reasoning, as do the 
snapshots presented in Fig. [T] Notice that in this region the membrane undergoes near- 
critical fluctuations. On the other hand, in the coexistence region, equilibrated membranes 
show complete phase separation and form large-scale lipid domains (Fig. [I]). 

The above analysis of the structure factors Si,(k) and Ss(k) and equilibrated membrane 
snapshots shows that in the part of the phase diagram where the lipid demixing curve closely 
approximates the binodal, the phase transition from the fluid state to the fluid-gel coexis- 
tence has a quasi- abrupt character. On the other hand, in the region where the lipid demixing 
curve strongly deviates from the binodal, the membrane shows characteristic near-critical 
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fluctuations (Fig. [I]), and the transition becomes continuous as the system approaches the 
critical point (DMPC/DSPC « 20:80, T c = 320.5 ±0.2 K - see [40J for details). Remarkably, 
this behavior is qualitatively similar to recent experimental observations on three-component 
lipid membranes [TT]: depending on the membrane composition, the transition to the two- 
phase coexistence takes place either in an abrupt manner - when the membrane does not pass 
through a critical point, or via critical fluctuations - when the membrane does pass through 
a critical point. A more detailed consideration of the phase diagram of the DMPC/DSPC 
system based on MC simulations will be published elsewhere [40J. 

In the fluid-gel phase coexistence region, large-scale phase separation takes place. Since 
no explicit or implicit penalties are imposed [4T1 - H3] on the domain size in our model, and 
domain growth is driven by minimization of the line tension energy, phase separation results 
in formation of a single circular-shaped domain of the minority phase (Fig. [I]). 

Effect of the cytoskeleton on phase separation in the membrane 

In the region of near-critical fluctuations the membrane is expected to be very sensitive to 
external perturbations, including the interaction with the membrane skeleton. It appears 
that in the vicinity of the critical point the interaction with the membrane skeleton leads to 
immediate condensation of the gel phase on the skeleton filaments, and formation of mem- 
brane domains, in a good agreement with results of experiments on lipid bilayers with a 
reconstituted actin skeleton [9]. These domains dynamically change their shape, but nev- 
ertheless stay pinned to the filaments (Fig. [2]). Notice that the effect of the filaments is 
remarkably robust with respect to the filament pinning density. We strongly believe that 
these domains represent the minimal model of membrane rafts [6j . The minimum character 
of this model stems from the fact that in the present scenario domain formation is thermo- 
dynamically driven, and does not require any active processes like chemical cross-linking of 
membrane components or lipid recycling. 

If the membrane is abruptly cooled down from the all-fluid state to a temperature in the 
fluid-gel coexistence region of the phase diagram, phase separation takes place, and domains 
of the fluid and gel phase start to nucleate and coarsen with time. In a free membrane, 
domains grow according to the power law R(t) ~ t n with the growth exponent n depending 
on the particular domain growth mechanism [44] . In our simulations for a free membrane, 
depending on the lipid composition and temperature, n takes values from 1/4 to 1/3 (Fig. 
[3]), consistent with general expectations for the domain growth in 2D [44] and in agreement 
with experimental results (see, e.g., [45]) (for a more detailed discussion, see [40]). 

In a finite-size system with a linear dimension L, domain coarsening eventually stops, 
and a single circular-shaped domain is produced (Fig. [I]) with the radius R^ ~ (X/7r) 1 / 2 L, 
where 0<X<l/2is the fraction of the minority phase. 

Remarkably, we found that the presence of the membrane skeleton strongly inhibits or 
even eventually prevents large-scale phase separation in the phase coexistence region (Fig. 
[2]). As a result, the radius of the membrane domains in this case is largely determined by the 
characteristic compartment radius R comp — £/2 of the filament network. The way we account 
for interaction of lipid molecules with filaments at their pinning sites is similar to a theoretical 
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Figure 2: Effect of the membrane skeleton on the phase separation in a DMPC/DSPC 
membrane. Representative snapshots of membrane configurations are shown for the free 
membrane (first row) and membrane interacting with a network of filaments at 100% (sec- 
ond row), 50% (third row) and 25% pinning density (fourth row). Snapshots for the free 
membrane, as well as for the membrane interacting with filaments at T = 321 and 322 K, 
represent fully equilibrated configurations; snapshots for the membrane interacting with fil- 
aments at T = 310 and 317.5 K correspond to equilibration time of 6 x 10 6 MC cycles (see 
text for discussion). For presentation purposes, the filaments are drawn thicker than they 
are in reality. Lattice size: 400 x 400; scale bar: 125 lattice units w 100 nm. 



model (the random-field Ising model [46]) implying the presence of static random position- 
dependent perturbations in a 2D system of spins. This model predicts that the initial power- 
law domain growth R(t) ~ t n is strongly slowed down at intermediate stages, and crosses over 
to an extremely slow logarithmic growth R(t) ~ logt; eventually the domains are expected 
to reach the perturbation strength-dependent maximum size R' < R^ in an exponential 
time to create an equilibrium disordered state [471 SH] . In our system, the perturbation 
strength is determined by the filament pinning density and the average compartment radius 
i? comp . Therefore, if the interaction of the membrane with the cytoskeleton is strong enough, 
and Rqq 2 3> i?c 0mp , one can expect that the domain growth stops when domains reach the 
characteristic size R coinp < R r < Roc- 

We indeed observed an extreme slowing down of the domain growth in the presence of 
the membrane skeleton and cross-over to the slow logarithmic growth (Fig. [3]). As is evident 
from the Figure, at the end of a simulation run of 6 x 10 6 MC cycles, the domain sizes are 
indeed R(t) ~ R comp - Although some growth is still observed at this stage, it is so slow, 
that at present it is unclear, whether the system will evolve toward a disordered equilibrium 
state featuring a number of small domains, as suggested in [471 SB], or the slow logarithmic 
growth will continue further to produce eventually (after an extremely long time) a single 
domain with the size R^. What is clear, though, that the time required for complete 
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Figure 3: Effect of the membrane skeleton on the the domain growth in DMPC/DSPC 50:50 
membrane abruptly cooled from the all-fluid state down to T = 310 K in the fluid-gel phase 
coexistence region. Kinetics of domain growth when the membrane is free (upper curve) 
and in the presence of membrane skeleton with the filaments pinning density of 50% (lower 
curve). Black solid line shows the power law dependence R(t) ~ t n with n = 0.32; dashed 
line shows the stage of the slow logarithmic growth R(t) ~ hit in the presence of membrane 
skeleton. Representative membrane configurations obtained in our MC simulations at time 
instants 10 4 , 10 5 , and 10 6 MC cycles are shown at the corresponding curves. Lattice size: 
400 x 400. 



equilibration is so long, that from the practical viewpoint one can state that indeed the 
presence of the membrane skeleton prevents large-scale phase separation in the membrane. 
This becomes especially clear if one takes into account that a membrane in a live cell is not 
in the equilibrium state, and a number of other processes affecting the membrane state and 
composition, e.g., lipid recycling, take place in parallel on different spatial and time scales. 

Thus, the cytoskelet on-induced inhibition of large-scale phase separation can serve as one 
of the possible explanations why micrometer-scale membrane domains are observed in giant 
plasma membrane vesicles [49], but not in cell membranes. 

In addition, our observations suggest that the established cryoprotective role played by 
the cytoskeleton [50] may consist in delaying phase separation-induced cold shock damage 
of living cells. 

Diffusion of lipids in the free membrane 

At higher temperatures, away from the phase transition and the near-critical fluctuations 
region, the membrane is in the homogeneous all-fluid state, and, not unexpectedly, lipid 
diffusion is normal (Fig. [4]). 

This picture changes upon approaching the critical point in the region of near-critical 
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Figure 4: Effect of the proximity to the phase transition on diffusion of DMPC lipids 
in a DMPC/DSPC 20:80 membrane. Mean-square displacement MSD(r) (a), local expo- 
nent /3msd(V) (b), time-dependent diffusion coefficient D(r) (c), and FCS autocorrelation 
G(r)/G(0) (d). Panels (a) - (c) top to bottom: T = 328, 322, and 321 K. Panel (d): 
T = 328 (left) and 321 K (right). For clarity, data for 322 K are omitted in panel (d). Panel 
(d) additionally shows fits to the FCS diffusion model Eq. Q giving /3 F cs = 1-01 at 328 K 

and /3pcs = 0.86 at 321 K ( ). For comparison, a fit of 321 K data with fixed /3pcs = 1.0 

is shown ( ). 



fluctuations where interpenetrating fluctuating domains form in the membrane. The anal- 
ysis of the immediate environment of the DMPC and DSPC lipids shows that, while the 
DSPC lipid does not have a pronounced preference for the state of its local environment 
and therefore is evenly distributed between fluid and gel domains, the DMPC lipid shows 
a strong preference for the fluid local environment and is thus predominantly partitioned 
into fluid domains. Under these conditions, the DSPC lipid shows no significant deviation 
from the normal diffusion behavior (data not shown). By contrast, quite unexpectedly, the 
DMPC lipid was found to demonstrate very pronounced anomalous subdiffusion (Fig. 

At a first glance, this subdiffusive behavior is rather surprising, especially in the light 
of the seminal work by Kawasaki [5]] demonstrating that the diffusion coefficient of a two- 
component mixture does not vanish upon approaching the critical point, which implies nor- 
mal diffusion, at least in the long-time asymptotic regime. A closer look at the data shows, 
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however, that the subdiffusive behavior is in fact transient, though it covers several orders 
of magnitude in time, and at longer times the cross-over to the normal diffusion takes place. 
One can also notice that at very short time lags, diffusion is normal as well. This behavior 
is also evident from the plots of the local exponent /3msd(^) and effective time-dependent 
diffusion coefficient D(r) extracted from the MSD{r) dependences (Fig. [4)3, c). 

How this cross-over from normal diffusion to a subdiffusive behavior and back to normal 
diffusion at long times can be explained? By recalling that the DMPC lipid is preferentially 
partitioned into fluid domains and having a look at the corresponding membrane snapshots 
(Fig. [T]), we realize that DMPC molecules diffuse on a dynamically rearranging fractal- 
like fluid phase pattern. Therefore, at early times, when the DMPC molecule explores its 
immediate environment, normal diffusion should take place. Later, though at times shorter 
than a characteristic time of rearrangement of this dynamic fractal structure, the motion 
of DMPC lipid molecules is subdiffusive, similar to what is expected in case of diffusion 
on fractal-like percolation clusters [52]. At much longer times, evolution of these clusters 
leads to dynamic percolation behavior (see, e.g., [53]). which results in a crossover from the 
subdiffusive motion to normal diffusion. In contrast to diffusion on static percolation clusters, 
where cross-over to normal diffusion occurs only above the percolation transition [52], in 
the case of lipid diffusion in a near-critical membrane, cross-over to normal diffusion will 
always take place irrespectively of whether the fluid phase is above or below the percolation 
threshold. 

This subdiffusive behavior is also clearly observed in our simulations of FCS experiments 
(Fig. [4]d) . Upon a closer approach to the critical point of the system, the autocorrelation 
functions of fluorescence intensity fluctuations G{r) progressively stronger deviate from the 
dependence expected for normal diffusion (i.e., Eq. Q with /3 F cs = 1). A good description 
of G(t) data in this case can only be obtained with /5pcs < 1- For example, we find that for 
DMPC/DSPC 20:80 mixture at T = 321.0 K the best fit is obtained with /3 FCS = 0.86 (Fig. 

and a smaller /3 F cs = 0.79 is required to fit G{r) at T = 320.7 K, which is closer to the 
critical point (data not shown). 

A new important result of the present work is that transient anomalous subdiffusion 
spanning several orders of magnitude in time can be observed in a system close to its critical 
point. We emphasize that the appearance of the transient subdiffusive behavior in the region 
of near-critical fluctuations is not related to any specific properties of the system and should 
be observed close to criticality in various systems independent of their origin, including, of 
course, multicomponent lipid membranes. 

It is well known that phase separation in lipid bilayers produces domains with typical sizes 
ranging from a few to several tens of micrometers (provided that there are no restrictions 
on the domain growth). In agreement with that, we observe in our simulations that within 
the coexistence region an equilibrated membrane always shows complete phase separation 
resulting in a single circular-shaped domain of the minority phase. Therefore, in this case, 
the only reasonable way to carry out FCS measurements requires parking the detection spot 
into the bulk of the majority phase away from the inter-phase boundary, or into the center 
of the minority phase domain. Exactly this approach is used in experimental FCS studies 
on membranes showing large-scale phase separation (see, e.g., [MH56] ). Not unexpectedly, 
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this approach results in normal diffusion in both phases with phase-dependent diffusion 
coefficients. This behavior is also observed in our simulations (data not shown). 

In contrast to that and quite surprisingly, an MC simulation study of Hac et al. [21] 
reported FCS curves strongly deviating from the normal diffusion model exactly in the 
phase coexistence region, i.e., where large-scale phase separation takes place. These results 
were attributed by the authors of [21] as resulting from subdiffusive motion of lipids, which 
is in contradiction with the macroscopic phase separation in the system. This seeming 
contradiction is resolved upon careful examination of the approach to simulations in [21]. 
There, the diameter of the FCS detection spot was approximately equal to the simulation box 
size. In this case, the FCS detection spot covers effectively the whole simulated membrane, 
and, even in the presence of large-scale phase separation, motion of lipids in both fluid and 
gel phases contributes to the FCS results, which explains deviations of FCS curves from 
the normal diffusion model observed in [21] in the fluid-gel phase coexistence region and 
rule out the interpretation of these results in terms of subdiffusion. Moreover, the situation 
simulated in [21] is very unlikely for experiments on single lipid bilayers within the phase 
coexistence region, where membrane domains typically have radii of a few to several tens of 
micrometers, whereas the typical FCS detection spot size is ~ 200 nm for the standard FCS 
[35], and down to 25 nm for STED FCS [36J. 

Here we should also mention a previous lattice-based simulation study [57] addressing the 
character of lateral diffusion of lipid molecules in DMPC/DSPC 50:50 membranes. There, in 
agreement with our results, diffusion is normal at T > 320 K, i.e. in the all-fluid membrane 
state (cf. Fig. [I]). Deviations from the normal diffusion for the DMPC/DSPC 50:50 mixture 
found in [57] within the temperature range T = 300 — 320 K (i.e. exactly where the large- 
scale phase separation for this lipid composition takes place - see Fig. [T]) were interpreted 
in terms of subdiffusive motion of lipids. As in [21], the spatial scales on which deviations 
from normal diffusion were observed in [57] approach the system size. Hopping of lipid 
molecules between the macroscopic gel and fluid domains definitely leads to deviations from 
normal diffusion at these large observation scales, although these deviations can hardly be 
interpreted in terms of anomalous subdiffusion. 

We emphasize that our results presented in this paper unambiguously show that local 
component and phase fluctuations in a near-critical membrane lead to the transient anoma- 
lous subdiffusion spanning several orders of magnitude in time. Under suitable experimental 
conditions, we believe, this subdiffusive behavior can be observed experimentally using, e.g., 
the (STED) FCS technique. 

Effects of membrane— cytoskeleton interaction on lipid diffusion 

Interaction with the membrane skeleton strongly affects the character of diffusion in the re- 
gion of near-critical fluctuations (Fig. [5]). With increasing the filament pinning density, the 
interaction of the membrane with filaments becomes more pronounced: it first enhances the 
anomalous diffusion and eventually leads to hop-diffusion of lipid molecules. In the presence 
of filaments, the faster diffusion process clearly corresponds to Brownian motion within com- 
partments defined by the membrane skeleton, whereas the slower diffusion is due to hopping 
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Figure 5: Effect of the membrane skeleton on diffusion of DMPC lipids in a DMPC/DSPC 
20:80 membrane at 322 K. Mean-square displacement MSD(r) (a), local exponent /3msd(V) 
(b), time-dependent diffusion coefficient D(r) (c), and FCS autocorrelation G(r)/G(0) (d). 
Panels (a) - (c) top to bottom, and panel (d) left to right: filament pinning density equals 
0% (free membrane), 25%, 50%, and 100%. 



between the compartments, and therefore it strongly depends on the filament pinning den- 
sity. Notice that gel phase condensation at the membrane skeleton substantially increases 
the effective thickness of filaments and thus enhances their influence on lipid diffusion. 

Remarkably, in this case our results for the mean-square displacement MSD{r) (Fig. 
[5^i) and effective time-dependent diffusion coefficient D{r) (Fig. [5^) are in good qualitative 
agreement with results of single-molecule tracking experiments in cell membranes [58] , where 
results were interpreted in terms of a picket-fence model. Our MSD(r) data also qualita- 
tively agree with a theoretical model for diffusion on an infinite periodic square meshwork 
of penetrable barriers [59] , except for in our case the short-time motion of a molecule within 
a compartment is subdiffusive due to near-critical fluctuations. We again emphasize that 
even a relatively low pinning density of the cytoskeleton filaments leads to appearance of 
rather strong barriers for diffusing lipids due to filament-induced local condensation of the 
gel phase. 

Our simulations of STED FCS experiments with the experimentally achievable detection 
spot size [361 EQ] show (Fig. [5]d) that, while at early times all FCS curves exhibit behavior 
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characteristic of subdiffusion, interaction with the membrane skeleton leads to appearance of 
a second slow diffusion component. The contribution of this slow component increases with 
the filament pinning density. Clearly, this component is related to the large-scale diffusion 
that involves crossing the filament-induced barriers on the membrane. This two-component 
character of G(t) is in qualitative agreement with the recent STED FCS measurements of 
lipid diffusion in cell membranes [60j . It would be very interesting to follow the dependence 
of the FCS autocorrelation functions and the respective apparent diffusion coefficients on 
the spatial scale using the variable detection spot FCS technique (HUES], as has been done 
in [60]. This work is currently in progress. 

Here, we point out once again that in our model the pinning sites do not represent 
obstacles for diffusing lipids, and affect lipid diffusion only indirectly, via the lower lipid 
mobility in the local gel-state environment (see Materials and Methods). We anticipate that 
in case where pinning sites would be presented by immobile particles preferentially wetted 
by gel-phase lipids, the effects of the membrane-filament interaction on diffusion of lipids 
should become even more pronounced. This issue will be addressed in our future work. 

Conclusion 

In conclusion, in this paper, we carried out Monte Carlo simulations of model two-component 
lipid membranes on experimentally relevant spatial scales and time intervals. This allowed 
us to better understand the details of the dynamics of phase and component separation, 
and by this means to address the relationship between membrane microheterogeneity and 
anomalous subdiffusion in cell membranes. 

We observed that interaction of the near-critical membrane with the cytoskeleton strongly 
affects phase separation: for the membrane in the state of near-critical fluctuations it leads 
to local phase separation and formation of dynamic domains with a size of a few tens of 
nanometers, which conform to the current understanding of membrane rafts [6]. 

We find that in a membrane showing near-critical fluctuations lipids show subdiffusive 
behavior covering several orders of magnitude in time. The interaction of the membrane 
with the cytoskeleton enhances subdiffusion and eventually leads to hop-diffusion of lipids. 

In the fluid-gel phase coexistence region of the phase diagram, interaction of the mem- 
brane with the cytoskeleton is found to tremendously slow down large scale phase separation 
in the membrane, which may explain the established connection between the cytoskeleton 
and cold shock resistance of organisms [50j . 

The concepts of lipid rafts [4J and cytoskeleton-related picket fence [5] are frequently dis- 
cussed as the alternative viewpoints on the origin of anomalous diffusion in cell membranes. 
In our paper, we bring these two concepts together to show that not only they do not contra- 
dict one another, but rather work in synergy, resulting in formation of cytoskeleton-induced 
dynamic lipid domains in the near-critical membrane. By this means we construct what, we 
believe, is a minimum raft model, since the domain formation is driven solely by thermody- 
namics and does not require either chemical cross-linking of membrane components, or lipid 
recycling. 
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Empirical phase diagram based on differential scanning calorimetry 
data 



A single-component membrane undergoes a phase transition from the gel state to the fluid 
state at its melting temperature. On the other hand, for a two-component lipid membrane 
the transition from the all-gel state to the all-fluid state is not immediate. (In fact, the 
two-component membrane undergoes two phase transitions: one from the gel state to the 
fluid-gel coexistence and, at a higher temperature, another one from the fluid-gel coexistence 
to the fluid state.) Therefore, from the practical point of view, one can speak about a broad- 
ened gel-fluid transition in a two-lipid system (compared to a single-lipid system) [ST14S3] . 
This broadened transition can be characterized by the onset and completion temperatures, 
which can be determined experimentally. By plotting the experimentally determined onset 
and completion temperatures as a function of the membrane composition, one obtains an 
empirical experimental phase diagram of the system. 

To determine the onset and completion temperatures from the experimental DSC data, 
the empirical tangent method was used. The outer slopes of the C(T) profiles of a range 
of compositions (DMPC/DSPC = 0/100, 10/90, 20/80, 90/10, 100/0) were fit with 
straight lines passing through the corresponding inflection points. The onset and completion 
temperatures in this case are defined as the respective intersections of the tangential lines 
with the zero-line. 

The experimental phase diagram of the DMPC/DSPC lipid mixture obtained in the 
present study is in a good agreement with the ones published earlier by other groups 



(Fig. SI) 



We remark here that evaluation of the onset and completion temperatures is not always 
unambiguous. In particular, our reanalysis of the experimental C(T) data from |S8j for 
DMPC/DSPC 20/80 and 10/90 mixtures gives onset temperatures lower than those reported 



in [S8] and very close to the ones we obtained from our data (Fig. SI). 



Monte Carlo simulations 

The fundamental step of the MC simulation consists of two sub-steps. In the first sub-step, 
an attempt is made to change the state of a randomly chosen lipid (from gel to fluid or vice 
versa). The second sub-step is an attempt to exchange the positions of a randomly chosen 
pair of next neighbors on the lattice. 

For each sub-step the change in the Gibbs free energy 

AG =ANf (AE 1 - TASi) + AiV 2 F (AE 2 - TAS 2 ) + AN^w^ + ANg F w^ 

+ AN^wf 2 G + ANlfw™ + AN?*w™ + AA^ 2 GF , (S1) 

is calculated. Here, AEi and ASi are the changes in the internal energy and entropy of 
a molecule of lipid i when it switches its state from gel to fluid, and w™ n are the next- 
neighbor interaction parameters of lipids i and j being in states n and m, respectively. In 
the simulation, the attempts of the state change and next-neighbor exchange are accepted 
with probability p = 1 for AG < and p = exp (-AG/RT) for AG > 0. 
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Figure SI: Empirical phase diagram of DMPC/DSPC constructed from excess heat capacity 
curves measured experimentally on multilamellar vesicle suspensions and calculated from 
MC simulation data. Experimental data: present work (•), Ref. [Si] (□), Ref. [S2j (o), 
Ref. [S4j (o), Ref. [S5j (v), Ref. [S6] (a), Ref. (STJ (<), Refs. (SSI ESI (>), our reanalysis 
of experimental C{T) data for DMPC/DSPC 20/80 and 10/90 from [S8] (>), Ref. [S9j (x). 
Monte Carlo simulation data: present work (■), Refs. [S3, S8J (+), Ref. [S9] (*). 



For an L x L lattice, one MC cycle consists of a chain of L 2 elementary MC steps. After 
each cycle, the enthalpy of the lattice is calculated as follows: 



H =N*AE 1 + N*AE 2 + N^wf? + Nl 



GF_GF , AT GF GF 
22 ^22 



N™w™ + N™v% + Ng*w™ + N?*w™ (S2) 



The excess heat capacity C (T) is calculated from the variance of equilibrium fluctuations 
of the total lattice enthalpy H as follows: 

C(T) = {{H-{H)f)/RT\ (S3) 

The lipid interaction parameters were adjusted using the approach previously described 
by Sugar et al [S3, ISlOj : temperature-dependent excess heat capacity curves were obtained 
from MC simulations for a range of membrane compositions and compared with experimen- 
tally measured heat capacity curves, and the parameters w™ n were varied until a reasonable 
agreement with experimental DSC data was achieved. Lipid interaction parameters for our 



model are listed in Table SI, Since a simpler lattice representation of the lipid system (lipid 
molecules arranged on a square lattice) was used in our simulations compared to the previ- 
ous studies [S3l IS8l - IS10| (lipids represented as dimers of acyl chains arranged on a triangular 
lattice) , not surprisingly, the lipid interaction parameters w^ n in our study differ from those 
used in previous publications. 



S3 



Table SI: Thermodynamic parameters of lattice-based MC simulations of the DMPC/DSPC 
lipid membranes used in the present work. 



Thermodynamic Parameters 


AE 1 

(Jmol 4 ) 


(Jmor 1 ) 


(Jmol^K" 1 ) 


AS 2 

(Jmol^K -1 ) 


(Jmor 1 ) 


(Jmol" 1 ) 


(Jmol" 1 ) 


(Jmol" 1 ) 


(Jmol -1 ) 


FF 
<2 

(Jmol" 1 ) 


26330 


50740 


88.653 


154.695 


1827 


1622 


4025 


4460 


1412 


502 



As is evident from Fig. SI , with this set of lipid interaction parameters, the heat capacity- 
based empirical phase diagram obtained in our MC simulations agrees well with experimental 
data (both ours and previously published), as well as with previously published results of 
lattice-based MC simulations of the same lipid system. 

For every lipid composition and temperature studied, the membrane was first equili- 
brated, if required. To do that, a system being initially at T = oc was brought to the 
equilibrium at a target temperature T by exercising the Monte Carlo procedure which ad- 
ditionally included a third MC substep consisting in position exchange of two randomly 
chosen lipid molecules. This approach, which was suggested in |Sllj and successfully ap- 
plied in [S3l IS8l IST2] is known to be very efficient in driving the system toward the equilibrium 
configuration. This procedure was propagated typically for 1.5 x 10 5 MC cycles, which is 
substantially longer than the typical time required by the total lattice enthalpy H to reach its 
equilibrium value at a given temperature in the presence of the random lipid exchange sub- 
step. After the completion of this procedure, the random lipid exchange substep is switched 
off, and the equilibrated system is ready for studies of its equilibrium properties. 

The simulation code was written in Fortran and compiled with Compaq Visual Fortran 
Compiler Version 6.6A (Compaq Computer Corporation, Houston, TX). Results were ana- 
lyzed using original dedicated routines written in MATLAB (The MathWorks, Nattick, MA). 
A reliable random number generator is essential for the success of an MC simulation. There- 
fore, we used the Mersenne Twister routine [S13J providing sequences of pseudo-random 
numbers equidistributed in 623 dimensions and characterized by an extremely long period of 
219937 _ ^ ~ 4 3 x io 6001 . Monte Carlo simulations were carried out on a workstation (Intel 
Core2 Quad Extreme CPU X9770 3.2 GHz, 4 GB RAM) running under Windows XP. Under 
these conditions, a simulation run on a 600 x 600 square lattice for 2 x 10 7 MC cycles takes 
about 700 h of CPU time. 



Calculation of static structure factors for lipids and lipid states 

To gain information on the character of spatial fluctuation of the local membrane composition 
and lipid state in the equilibrium, the static structure factors for lipids (DMPC and DSPC) 
Si,(k) and lipid states (fluid and gel) Ss(k) were calculated independently. 
The structure factor 

S(k) = l + pg(k) (S4) 



S4 



is expressed via the Fourier transform g(k) of the pair distribution function of particles 



where p is the spatially averaged number density of particles |S14j . To avoid artifacts in 
spatial and angular averaging, the fact that the particles cannot take arbitrary positions, 
but rather occupy sites on a square lattice, was taken into account. 
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